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We set up a framework in which in-medium charmonium properties are constrained by thermal 
lattice QCD and subsequently implemented into a thermal rate equation enabling the comparison 
with experimental data in heavy-ion collisions. Specifically, we evaluate phenomenological conse- 
quences for charmonium production originating from two different scenarios in which either the 
free or the internal energy are identified with the in-medium 2-body potential between charm and 
anti-charm quarks. These two scenarios represent J/ip "melting temperatures" of approximately 
1.25 T c ("weak binding") and 2T C ("strong binding"), respectively. Within current uncertainties in 
dissociation rates and charm-quark momentum spectra, both scenarios can reproduce the centrality 
dependence of inclusive J/tp yields in nuclear collisions at SPS and RHIC reasonably well. However, 
the "strong-binding" scenario associated the the internal energy as the potential tends to better 
reproduce current data on transverse momentum spectra at both SPS and RHIC. 



PACS numbers: 25.75.-q, 12.38.Mh, 14.40.Lb 



I. INTRODUCTION 

Systematic studies of charmonium production in ul- 
trarelativistic heavy-ion collisions (URHICs) are hoped 
to reveal the properties charm-anticharm quark bound 
states in hot and dense matter. A central goal is to utilize 
these insights in the search for the Quark-Gluon Plasma 
(QGP) and infer some of its basic properties, such as 
color screening, flavor transport and the relevant degrees 
of freedom (see, e.g., Refs. pHS] for recent reviews). To 
carry out this program quantitatively requires the com- 
bination of theoretical approaches to evaluate charmo- 
nium spectral function in medium with phenomenolog- 
ical models furnishing a realistic context for computing 
observables. The former category includes first-principle 
lattice-QCD (1QCD) calculations of charmonium 

correlation functions and heavy-quark (HQ) free ener- 
gies, as well the recently revived potential models [8hT5| 
to compute the charmonium spectrum in the QGP. In 
the latter category, kinetic rate equations, transport and 
statistical models have been pursued [i~6H23l ]. However, 
there are currently rather few calculations with quantita- 
tive connections between these two categories (see, e.g., 
Ref . [24[ ) ■ It is the purpose of the present paper to elab- 
orate such connections. 

The phenomenological building block in our work con- 
sists of a kinetic rate equation which accounts for char- 
monium dissociation and regeneration in a thermally ex- 
panding hreball. The in-medium charmonium properties 
figuring into the rate equation are inferred from spectral 
functions which are constrained by euclidean correlators 
and HQ free energies computed in thermal 1QCD. The 
spectral functions are adopted from a recent potential- 
model calculation using a thermodynamic T-matrix [TEj 
with HQ free (or internal) energies as in-medium driv- 
ing kernel. From these spectral functions we extract in- 
medium binding energies and bound-state masses which 
determine the reaction rates and equilibrium limit for 



the rate equation. We explicitly verify that the "re- 
constructed" spectral functions yield euclidean correla- 
tor ratios which vary little with temperature (say, within 
±10%), as found in 1QCD. In view of the ongoing de- 
bate as to whether the free or internal energy is a more 
suitable quantity to be identified with a HQ potential, 
we will explore both possibilities as supplied through 
the T-matrix calculations in Ref. As found there 

and in previous works the free and internal en- 

ergy may be considered as providing an lower and up- 
per limit, respectively, on the dissociation temperatures 
of charmonia, which in the following we will refer to as 
weak- and strong-binding scenarios. We will elaborate 
the phenomenological consequences of both scenarios for 
charmonium observables in URHICs at the Super Proton 
Synchrotron (SPS) and at the Relativistic Heavy-Ion Col- 
lider (RHIC), including centrality, transverse- momentum 
(pt) and rapidity (y) dependencies. 

This paper is organized as follows: In Sec. [H] we ex- 
amine the equilibrium properties of charmonia obtained 
from 1QCD via the potential model and extract the key 
parameters needed for the kinetic approach of charmo- 
nium production in URHICs. In Sec. IIIII we review the 
required ingredients of our previously constructed rate 
equation for calculating the inclusive yield as well as the 
transverse-momentum spectra of charmonia. In Sec. IIVI 
we compare our numerical results with experimental data 
at SPS and RHIC. We conclude and give an outlook in 
SecM 



II. EQUILIBRIUM PROPERTIES OF 
CHARMONIA 



In this section we discuss the equilibrium properties 
of charmonia (^=J/ip, Xc and tp') in the QGP medium 
as calculated from 1QCD and extract the temperature- 
dependent quantities needed in the kinetic approach. 
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Specifically, these are the in-medium charm-quark and 
charmonium masses, m* and m^, respectively, the char- 
monium dissociation rate, r^ lss , and the charmonium dis- 
sociation temperature, T^ iss . The latter defines the tem- 
perature above which regeneration is not operative (i.e., 
the gain term of the rate equation is set to zero). We 
define the binding energy as 



e B (T) = 2m*(T)-m«(T) 



(1) 



In Sec. Ill Al we review basic relations between correla- 
tors and spectral functions and extract pertinent quanti- 
ties from the T-matrix approach of Ref. [Hj]. In Sec. Ill Bl 
we detail our calculations of the charmonium dissocia- 
tion widths which we assume to equal the total width of 
the charmonium. In Sec. Ill CI we combine the extracted 
quantities to "re-construct" charmonium spectral func- 
tions and evaluate corresponding euclidean correlator ra- 
tios in light of 1QCD results. 



A. Euclidean Correlators and Potential Model 

Lattice QCD currently provides equilibrium properties 
of charmonium mainly in terms of two quantities: 
1. The free energy, Fqq(t;T), of a static pair of heavy 
quark and anti-quark. This is the main input for re- 
cent potential models. It remains controversial to date 
whether the free energy, the internal energy, 



U, 



Q Q (r;T) = F Q 



}(r;T)-T- 



9F Q o(r;T) 



dT 



(2) 



or any combination thereof, should be identified with a 
static QQ potential at finite temperature, T. In this 
work we therefore study two scenarios where either Uqq 
or Fqq is used as potential, Vqq. Since the internal 
energy leads to stronger binding than the free energy, we 
refer to the former and latter as strong- and weak-binding 
scenario, respectively. 

2. The two-point correlation function of a quarkonium 
current, j a , with hadronic quantum number a, 



G Q (r,r)HO Q (r,r)jt(0,0))) , 



(3) 



computed as a function of imaginary (euclidean) time, r 
(also called temporal correlator). The imaginary part 
of the Fourier transform of the correlation function, 
G a (r, r) , is commonly referred to as the spectral func- 
tion, 



cr a (uj,p) = --ImG a (w,p) , 
which is related to the temporal correlator via 



(4) 



G a (r,T)= / dija a {u,T)K(Lj,T,T) (5) 
Jo 



with the finitc-T kernel 



K{u,t,T) 



cosh[(w(r - 1/2T)] 
sinh[u;/2T] 



(6) 



Lattice QCD results for two-point correlation functions 
are usually normalized to a "reconstructed" correlator 
evaluated with the kernel at temperature T, 

/>oo 

GT(r,T) = / du a a (u,T*) K(w,r,T) , (7) 
Jo 

but with a spectral function at low temperature, T*, 
where no significant medium effects are expected. The 
correlator ratio, 



R a (r,T) = G a (r,T)/GT(r,T) , 



(8) 



is then an indicator of medium effects in G a (r, T) 
through deviations from one. Current 1QCD calcula- 
tions find that the correlator ratio, R a (r, T), in the pseu- 
doscalar (rj c ) and vector (J/i/j) channel are close to 1 
(within ca. 10%) at temperatures up to 2-3 T c pi ua. |26| . 
In the P-wave channels (scalar and axialvector) the cor- 
relators ratios are substantially enhanced over 1 at large 
r. This feature is believed to be due to "zero-mode" con- 
tributions (at uj—0) which are related to the scattering 
of a charm (or anti-charm) quark, c — > c (or c — > c), 
rather than to cc bound-state properties [27[ ■ This inter- 
pretation is supported by studies of the T-derivative of 
P-wave correlator ratios, which exhibits a much smaller 
variation (in the limit that the zero-mode part is a 8- 
function, <7 zm (w) oc 8(uj), its contribution to the temporal 
correlator is a constant) [28l.[29j. 

In principle, the in-medium properties of charmonia, 
such as pole mass, in-medium width and dissociation 
temperature, are fully encoded in their spectral func- 
tion. However, the finite number of data points for the 
two-point correlator computed in 1QCD severely ham- 
pers the inversion of the transform in Eq. ([5]), rendering 
the determination of the spectral function difficult. In 
fact, potential-model analyses have shown that the use 
of either the free or internal energy can lead to agree- 
ment with 1QCD correlators, albeit with rather differ- 
ent underlying binding properties (and associated dis- 
sociation temperatures). In the present work we there- 
fore adopt the following strategy: we first extract the 
charm-quark masses and charmonium pole mass from a 
potential model, allowing us to define the binding en- 
ergy according to Eq. ([1]). The latter is an important 
ingredient in the quantitative evaluation of the ^ disso- 
ciation rate (30j , which will be done in the following sec- 
tion pTB) in a perturbative ( "quasifree" ) approximation. 
We then "reconstruct" in-medium charmonium spectral 
functions using a relativistic Breit-Wigner + continuum 
ansatz, where the $ width and mass figure into the Breit- 
Wigner part while the continuum is determined by the 
open-charm threshold (2m*). For a more realistic eval- 
uation, we include polestrength factor, Zy(T), for the 
Breit-Wigner strength and a nonperturbative rescatter- 
ing enhancement in the continuum [ill . [l3[ . The vanish- 
ing of the polestrength factor furthermore serves to esti- 
mate the dissociation temperature of the ground state in 
each channel. 
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FIG. 1: (Color online) Temperature dependence of in-medium 
charm quark mass. The solid (dashed) line is for the strong 
(weak) binding scenario. 



For definiteness we employ the potential model of 
Ref. [l5l | where quarkonium spectral functions and corre- 
lators have been calculated in a thermodynamic T-matrix 
approach, consistent with vacuum spectroscopy and in- 
cluding relativistic corrections for a proper description of 
scattering states. The calculations in there have been 
carried out for both free and internal energies as po- 
tential, and for two different 1QCD inputs 0, H|. In 
both cases (and for both potentials), an approximate con- 
stancy (within ±15%) of the correlator ratios for pseu- 
doscalar charmonium has been found (see lower panels 
of Fig. 12 and 14 in Ref. [Hj]). We believe that these 
results provide a reasonable representation and bracket 
for potential-model results. In Fig. Q] the temperature- 
dependent charm-quark mass is displayed, which is iden- 
tified with the asymptotic value of the HQ potential, 



m *(T) = m c + V QQ (r^^-T)/2 . (9) 



The in-medium masses decrease with temperature ap- 
preciably, while the magnitude of m* (T) is significantly 
smaller in the weak-binding compared to the strong- 
binding scenario. As expected, the binding energies 
(plotted in Fig. [5]) also decrease with T, again being sig- 
nificantly smaller in the weak-binding scenario. These 
features are, in fact, the main reason that both scenar- 
ios can be compatible with the small variations found in 
the 1QCD correlator ratios: for weak/strong binding, a 
small/large constituent mass combines with a small/large 
binding energy, respectively, leading to an approximate 
compensation in the bound-state mass, to$(T), recall 
Eq.Q. 




T (GeV) 

FIG. 2: (Color online) Temperature dependence of J/tj) bind- 
ing energy. The solid (dashed) line is for the strong (weak) 
binding scenario with V q q=U q q{F q q). 



B. Dissociation Rates 

The inelastic dissociation rate of charmonia in the 
medium, r^ lss (T), plays a central role in URHIC phe- 
nomenology as it directly governs the time dependence 
of their abundance in an underlying rate equation (see, 
e.g., Eq. (|2"3")l below). It contributes to the total width of 
the corresponding charmonium spectral function (in ad- 
dition to elastic scattering which we neglect in the present 
work). However, it is currently not possible to quantita- 
tively extract the width from temporal correlators f3lj ]: 
for phenomenologically relevant values of the width of a 
few tens of MeV (or even up to 200 MeV) the correla- 
tor ratios are affected at a level of ~5% [ll|, which is 
too small to be discerned from other uncertainties at this 
point. Furthermore, elastic and inelastic collisions con- 
tribute to the width with an a priori unknown partition 
(for relatively loosely bound states one expects the dis- 
sociation width to be dominant). In the present work we 
therefore calculate the dissociation rate based on pertur- 
bative QCD with an effective strong coupling constant 
which we later fine tune to reproduce the observed J/ip 
suppression in central A-A collisions. The resulting value 
of a s turns out to be ~0.3, quite compatible with the 
short-distance (color-Coulomb) term in the effective po- 
tential used to extract the binding energies discussed in 
the previous section. 

Let us briefly discuss basic mechanisms for charmo- 
nium dissociation in medium. In the QGP, the leading- 
order process is naively given by the gluo-dissociation 
process introduced more than 30 years ago |32| , & + g — > 
c + c. 1 However, for small binding energies this process 



In the language of effective field theory (EFT) [1411 . this corre- 
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becomes inefficient (due to a shrinking phase space) and 
is superseded by inelastic reactions with an extra par- 
ton in the final state, i + ^—ti + c + c (i=g,q,q) |30||. 2 
Following Ref. [30j, we treat this process as a "quasi- 
elastic" collision between the parton i from the medium 
and the core quark in the bound state, i + d — > i + c. 
In this quasifree approximation the charmonium dissoci- 
ation cross section becomes twice the elastic cross section 
between the parton and the c quark, a^ ss = 2af i: as, e.g., 
given in Ref. [33j]. To account for the leading kinematic 
correction from the residual binding energy, the incom- 
ing parton needs to be energetic enough to break up the 
bound state, which sets a lower limit for the incoming 
parton momentum. Overall 4-momentum conservation 
for the process i + $ — > i + c + c is maintained by as- 
signing the binding energy to a decrease in mass of the 
initial-state charm-quark, d, i.e., m c i = m c — eb- In 
addition, we have introduced a Debye mass, ran = gT , 
into the denominator of t-channel gluon-exchange prop- 
agator, 1/t — > l/(t — m 2 D ), to regulate the divergence for 
forward scattering (the strong coupling in is taken 
consistently with the coupling constant a s ). 

In the hadron gas (HG) phase, we employ a flavor- 
SU(4) effective Lagrangian approach [34], [35[ to estimate 
the inelastic cross sections with tt and p mesons. As we 
will see below most of the charmonium dissociation in 
the hot medium occurs in the QGP, because the typical 
density of hadrons in HG is much smaller than that of 
partons in QGP, while the dissociation cross sections are 
of similar magnitude (around lmb). 

The dissociation rate of a charmonium state at finite 
3-momentum, p, can be obtained from the inelastic cross 
section by a convolution over the thermal distribution, 
f l (uj;T), of medium particles in the QGP or HG, 



rf s ( P ,T) = 



d 3 k 
J2nf 



( 10 ) 

Here, v Te \ is the relative velocity of ^ and medium par- 
ticle i. The temperature dependence of the quasifree 
dissociation rates for the J/ip and \c m the QGP are 
plotted at p—0 vs. T in the upper panel of Fig. [31 and 
for the J/ip as a function of p at selected temperatures 
in the lower panel. In the weak-binding scenario, there 
is rather little difference between the dissociation rates 
of J/ip and Xc especially above T = 200 MeV. Only in 
the strong-binding scenario the larger J/ip binding en- 
ergy makes a large difference, suppressing its destruction 



sponds to the color-singlet to -octet transition, albeit final-state 
interactions in the c + c octet state are neglected here. However, 
they are repulsive and of order 0(1/N^) which renders them 
numerically very small. 

In EFT language, some of these processes (e.g., the ones involv- 
ing t-channel gluon exchange between the thermal parton i and 
one of the charm quarks in the bound state) correspond to the 
Landau damping contribution to the width. 
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FIG. 3: (Color online) Upper panel: temperature dependence 
of dissociation rates for J /if) and \c calculated in quasifree 
approximation for the strong-binding (solid line: J/ip, dot- 
dashed line: Xc) and weak-binding scenarios (dashed line: 
J/ip, dash-double-dotted line: Xc)- The dotted line indicates 
the contribution from J/ip gluo-dissociation rate in the strong 
binding scenario using a s =0.32 in the underlying cross sec- 
tion expression [3^ ]. Lower panel: 3-momentum dependence 
of the J/ip dissociation rate in the strong- and weak-binding 
scenarios (solid and dashed lines, respectively). 



by, e.g., a factor of ^5 at T ~200MeV relative to the 
Xc and ip' (not shown); this difference becomes larger 
(smaller) at smaller (larger) T. For comparison we also 
calculated the rate due to the gluo-dissociation mecha- 
nism employing the expression derived in Ref. 32j with 
the same a s =0.32 as in the quasifree rate and with es 
obtained from the strong-binding scenario (note that the 
Coulombic binding is much smaller), which turns out 
to be inefficient for dissociating J/^'s (and even more 
so for the excited states) and is thus neglected in the 
following. The 3-momentum dependence of the rates 
shows a monotonous increase with increasing p, which 
becomes more pronounced with increasing binding en- 
ergy (for larger £j a finite 3-momentum facilitates the 
break-up since, on average, a larger center-of-mass en- 
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ergy is available in the collision of the bound state with 
thermal partons). This increase is a simple kinematic 
consequence of a monotonously increasing (or even con- 
stant) cross section with finite threshold and an increas- 
ing parton flux encountered by the moving J/ip. 

In the following section we determine the interplay of 
the charmonium widths with their binding energies ex- 
tracted in the previous section in order to assess their 
survival in the QGP (e.g., as a resonance). 



C. Spectral Functions 

We now turn to the construction of the charmonium 
spectral functions and the corresponding two-point cor- 
relation functions, and constrain the latter by 1QCD 
computations. In doing so we introduce an additional 
quantity, namely the polcstrcngth, Z^(T), of the Brcit- 
Wigner bound-state part characterizing the disappear- 
ance of the state from the spectrum. The vanishing of 
Zij,(T) will be used to characterize the dissolution tem- 
perature above which regeneration in the rate equation 
is inoperative. 

We first construct a model spectral function in vacuum, 
consisting of a zero- width bound-state and a perturbativc 
(leading order) continuum part, 



cr^(w) = Ay S(lu — mqj) 
B*N, 



+ 



8tt 2 



8(w-« )uVl 



2 



(a + b- 



(11) 



Here, N c —3 is the number of colors and the coefficients 
(a, b) = (1, —1), (2, 1) characterize the scalar and vector 
channel, respectively @. The open-charm threshold in 
vacuum, sq, is assumed to be given by twice the free D- 
meson mass, sq = 2mjj = 3.74 GeV. The coefficient A is 
related to the overlap of the wave- function, Rj/^(0), or 
its derivative, R' (0), at the origin [9l. [36|. 



A 



3N C 
2ir 



Ac 



(12) 

These quantities can be estimated from the electromag- 
netic decays widths via 36] 



r — 

-L P.P. 



ie 2 Q a 2 N c 



lUej,a 2 N c 



\R JN m 2 , r 77 = \ -|i?^(0)| 2 

(13) 

where a— 1/137 is the electromagnetic coupling constant 
and eQ = 2/3 the charge of the charm quark (we use 
r ee =5.55keV for the J/tp and T 77 =2.40keV for the Xco)- 
The resulting relations between Ay and r^_j. eei77 are 



81^ 2 //y, 

327TO 2 



r,. 



81^Lo r 
12SW 77 • 



(14) 



The J/ip and Xco masses are taken at their empirical vac- 
uum values. The coefficient By in the continuum part 
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FIG. 4: (Color online) Temperature dependence of the 
strength of the resonance part of the S-wave spectral func- 
tion, Z<s,(T). The solid (dashed) line is for the strong (weak) 
binding scenario. 



of Eq. (fTTj) equals one in the non-interacting limit. To 
account for rescattering, which is particularly important 
close to threshold, we scale it up to match the contin- 
uum as calculated from the vacuum T-matrix in Ref . [15| , 
amounting to Bj/q, ~ 2 and B Xc ~ 4 in the vector and 
scalar channel, respectively. For simplicity we neglect tp', 
X' c and higher excited states which play little role in the 
correlator ratios. 

At finite temperature we replace the <5-function bound- 
state part by a relativistic Breit-Wigner (RBW) distri- 
bution while the continuum part is assumed to be of the 
same form as in the vacuum, 



ay (00) = Ay Zy(T) 



2uo 



*y(T) 



7T (w 2 -m!(T)) 2 + W 2 r*(T) 2 



ByN c 
8tt 2 



G(w- s(T))uj\ 1 



S (T) 2 



(15) 



The in-medium continuum edge, s(T), is now taken as 
the charm-quark threshold at finite temperature, s(T) = 
2m* (T) , consistent with the potential model, see Fig. Q] 
The RBW term includes the in-medium charmonium 
mass, my(T), extracted from Eq. (JlJ based on Figs. Q] 
and [2] the width identified with the inelastic disso- 
ciation width discussed in the previous section, and the 
aforementioned polestrength factor, Zy(T). The latter is 
adjusted to minimize the deviation of the correlator ra- 
tios from one. The resulting Z(T) for J/tp (vector chan- 
nel) is plotted in Fig. 21 from which we extract its dissoci- 
ation temperature Tj¥^=2.0(1.25)T c in the strong (weak) 
binding scenario. Similar analysis in the scalar channel 
yields \c dissociation temperatures of T^ SS =1.3(1.0)T C 
in the strong (weak) binding scenarios. We assume that 
Xci and Xc2 have the same dissociation temperatures as 
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FIG. 5: (Color online) Spectral functions in the vector chan- 
nel. The upper (lower) panel is for the strong (weak) binding 
scenario. 



FIG. 6: (Color online) Ratio of vector channel correlator to 
the reconstructed correlator. The upper (lower) panel is for 
the strong (weak) binding scenario. 



the Xco- F° r "0' we simply assume its dissociation tem- 
perature to be T c for both the strong- and weak-binding 
scenarios. 

To comprehensively illustrate the medium effects we 
plot the final spectral functions for the vector channel 
in the strong- and weak-binding scenario in the QGP in 
Fig. O and their corresponding correlator ratios in Fig. El 
the spectral functions for the scalar channel are displayed 
in Fig. [7] We see that the correlator ratios are indeed 
close to one, as found in 1QCD. In the hadronic phase (not 
shown), we assume vacuum masses for both charmonia 
and open-charm hadrons, which automatically ensures 
that the correlator ratios are close to one (deviations due 
to small charmonium widths in hadronic matter are neg- 
ligible) . 

We are now in position to implement the in-medium 
properties of the charmonia into a kinetic rate equation 
in a thermal background, paving the way for applications 
to experimental data in heavy-ion collisions. 



III. KINETIC APPROACH 

The evolution of charmonium yields and spectra in a 
nucleus-nucleus {A- A) collision can be roughly divided 
into three stages. In the first, "primordial" stage (at 
small times r ~ 0) charm-quark pairs are produced in 
initial hard nucleon-nucleon (N-N) collisions. In the sec- 
ond, "pre-equilibrium" stage (0 < r < r ) charmonia 
are typically in the formation phase, but the so-called 
pre-resonance states are already subject to dissociation 
by passing- by nucleons. The third, "equilibrium" stage 
( r ^ T o) starts once the hot and dense medium has 
thermalized and lasts until thermal freezeout after which 
hadrons stream freely to the detectors. During this stage 
charmonia are subject to dissociation by particles of the 
hcatbath, but detailed balance requires that cc pairs in 
the medium can also recombine and "regenerate" char- 
monia (below their respective dissociation temperature) . 

Guided by the above scheme, we have organized the 
remainder of this section as follows: In Sec. IIII Al we dis- 
cuss primordial charmonium production and its modi- 
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FIG. 7: (Color online) Spectral functions in the scalar chan- 
nel. The upper (lower) panel is for the strong (weak) binding 
scenario. In the weak-binding scenario XcO has already melted 
at T c . 



fixations in the pre-equilibrium stage, including rapidity 
and transverse- momentum dependencies. In Sec. lIII Bl wc 
briefly recapitulate a simple fireball model which serves 
as a thermally evolving background for the charmonium 
rate equation. The latter is introduced in Sec. IIII CI 
along with its main ingredients related to in-medium 
open-charm and charmonium properties. In Sec. IIII Dl 
we specifically address the transverse-momentum de- 
pendence of charmonium production in the equilibrium 
stage. 



A. Cold Nuclear Matter Effects 

Charmonium production yields and pt distributions 
following from the primordial and pre-equilibrium stage 
form the initial conditions for the thermal rate equation 
discussed in the next section. The starting point are pri- 
mordial spectra based on experimental data in p-p colli- 
sions, modified by corrections specific to proton-nucleus 
(p-A) and extrapolated to A- A collisions. The corrections 



are commonly referred to as cold-nuclear-matter (CNM) 
effects, which include the following: 

1) Nuclear shadowing, the modification of the initial par- 
ton distribution functions in a nucleus relative to those 
in a proton, which affects both open and hidden charm 
yields in A- A relative to p-p collisions. 

2) Cronin effect, the increase of the mean p t of produced 
charmonia in A- A relative to p-p associated with initial- 
state parton scattering prior to the hard production. 

3) Nuclear absorption, the dissociation of prc- 
charmonium states by passing-by nucleons. 

In the present work, we assume a factorization of the 
charmonia phase-space distribution function into spatial 
and momentum parts, 



h(b,x t ,p t ,T Q ) = fe(b,x t ,Ta) fa(b,Pt,T ) 



(16) 



(6: impact parameter of the A- A collision, which we use 
as a measure of centrality equivalent to iV part , the num- 
ber of nucleon participants). The Cronin effect is read- 
ily implemented into the 3-momentum dependent part 
f<s{b,Pt,To) via a Gaussian smearing of the charmonium 
Pt distribution in p-p collisions, f^, p (pt), 



far(b,pt,T ) 



d 2 q t 



■ exp 



Qt 



-3*1) ■ 

The nuclear increase of the average p\ , (Apf) = {p\)aa — 
(Pt)ppi 1S estimated within a random- walk treatment of 
parton-nucleon collisions [37l | as being proportional to 
the mean parton path length, (l), in the cold medium: 
(Ap 2 ) = a g N (I). The coefficient a g N is estimated from 
p-A data at SPS || and d-Au data at RHIC [H . We 
use a sJV =0.076GeV 2 /fm for ,/I=17.3 AGeV Pb-Pb colli- 
sions and a gAr =0.1(0.2) GeV 2 /fm for ^=200 AGeV Au- 
Au collisions at mid (forward) rapidity. 

Nuclear absorption is evaluated within a Glauber 
model; the resulting charmonium distribution in the 
transverse plane takes the form 



f 9 (b,x t ,r ) = Ay- 



dy 



dz dz'pA (x t ,z) p B (b - f t , z') 



x exp I - I dz A p A {xt,z A )(T ahs 



x exp \ ~ J dz B p B {b - x t ,z B )a- ahB 

z' > 

(18) 

where pa,b are Woods-Saxon profiles of nuclei A and 
B and Ay=1.8 represents the rapidity coverage of our 
thermal fireball (see Sec. IIII Bl below). For the pp char- 
monium production cross per unit rapidity we take the 
values d<r* /dy=37nb for ^1=17.3 AGeV Pb-Pb (U 
(with ca. 40% uncertainty) and dcr^,/d?/=750(500) nb for 
yi=200AGeV Au-Au [42| at mid and forward rapidity 
(with ca. 10(20)% uncertainty). We utilize an effective 



8 



^-N absorption cross section, tx a b s , to parameterize both 
nuclear shadowing and absorption. Applying Eq. (fT5j) 

to p-A collisions at SPS we obtain u J J^=7.2>±\ mb from 
the recent NA60 data at £i a b=158GeV (corresponding 
to y/sNN~17.3 GeV) (43|. The new experimental mea- 
surement at 158 GeV turns out to give a significantly 
larger value than previously available for 400 GeV pro- 
ton projectiles, <7 a b s — 4.4mb [3] (the latter has been 
confirmed by NA60 14310, which has been used in our 
previous calculation s 12 ll 1 2 31 . The comparison with re- 
cent PHENIX data HaWyields cr abs ~ 3.5 mb (5.5 mb) 
for -^=200 AGeV Au-Au collisions at mid rapidity, \y\ < 
0.35 (forward rapidity, \y\ £ [1.2,2.2]). For simplicity, we 
assume the same absorption cross sections for the \c as 
for the J/ip. However, for excited states er a bs is expected 
to be significantly larger, even if they are not fully formed 
when the dissociation occurs. Taking guidance from 
the NA50 measurement with 400 GeV protons, we use 

crfbs - 13 mb at \/s=17.3 AGeV and o^ s ~ 6.5(10) mb 
at y / s=200AGeV for mid (forward) y. For each char- 
monium state, the number surviving the pre-equilibrium 
stage in an A-A collision at impact parameter b thus 
amounts to 



Ny (b) 



fy(b,Xt,T )d 2 x t . 



(19) 



The rather pronounced rapidity dependence of <7 a b s at 
RHIC casts doubt on interpreting this quantity as an 
actual absorption cross section. It seems more reason- 
able to associate its increase at forward y with nuclear 
shadowing |46j since the dissociation kinematics is very 
similar between mid and forward rapidity. While this 
does not affect the use of our "effective" cr a b s , it does im- 
ply a nuclear shadowing effect on the open-charm cross 
section in A-A collisions (which is an important ingredi- 
ent in the calculation of regeneration). As a "minimal" 
scheme we therefore associate the additional absorption 
of the J/ip yield at forward y (relative to mid rapidity) 
with a suppression of open charm production caused by 
shadowing, while we assume no shadowing corrections at 
mid rapidity. Thus, at both SPS and RHIC the num- 
ber of primordially produced cc pairs at mid-rapidity is 
calculated from the pp cross section as 



■ j dcr^ 
N™ d (b) = Ay^ 



Tab (b) 



while for forward y at RHIC we use 

da: 



N£(b) = Ay- 



)>ii 



Ay 



gfor 

TA B (b)-^ 

y=1.7 u nuc 



(20) 



(21) 



Here, Tab (b) is the usual nuclear overlap function 
and 5 nuc denotes the J/i/j suppression factor due to 
CNM effects, parameterized by cr a b s in the Glauber for- 
mula, Eq. ([Tg]). In particular, the ratio S^JS^ rep- 
resents the extra suppression associated with nuclear 
shadowing, operative for both J/ip and cc production. 



The input charm-quark cross section in pp is taken 
as ddcc/dy (j/=0)=2.2 fib at SPS (according to the re- 
cent compilation of data in Ref. |47|), and as da^c/dy 
(y=0)=123±40 ub at RHIC (in line with recent PHENIX 
measurements 48]). At forward rapidity we assume the 
pp charm-quark cross section to be reduced by 1/3 ac- 
cording to recent measurements (49| . 



B. Fireball Model 

Once the nuclear-collision system thermalizes, its tem- 
perature is the key quantity connecting to the in-medium 
properties of the charmonia as discussed in Sec. |TTJ As 
in our previous work [2lj. [23l |30| , we estimate the time 
evolution of the temperature using an isentropically and 
cylindrically expanding isotropic fireball characterized by 
an eigenvolume, 

Vfb(t) = (an, + v z t + \a z r 2 ) tt (Rq + \ai_T 2 ) 2 . (22) 

The fireball expansion parameters v z , a z , a ± are chosen 
such that the hadron spectra at thermal freezeout are 
consistent with the empirically extracted light-hadron 
flow in resemblance of hydrodynamical calculations (we 
actually use a relativistic form of a± (r) which limits the 
surface speed, v s (t)=ci±t, to below c). The initial trans- 
verse radius i?o represents the initial transverse overlap 
of the two colliding nuclei at a given impact parameter 
b, while the initial longitudinal length, zq, is related to 
thermalization time tq through zq ~ Az/to where Ay=1.8 
represents the typical longitudinal rapidity coverage of a 
thermal fireball. We assume that at a formation time 
of r = 1.0(0.6) fm/c the medium at SPS (RHIC) first 
thermalizes with all the entropy, <Stot(&), being built up. 
The latter is estimated from the multiplicities of observed 
charged particles and assumed to be conserved during the 
adiabatic expansion. We can then use the entropy den- 
sity, s(r) =Stot/^FB at each moment r in the evolution to 
infer the fireball temperature once we specify the equa- 
tion of state of the medium, specifically s(T). The QGP 
is modeled by an ideal gas of massive quarks and glu- 
ons while the hadronic phase is approximated by a non- 
interacting resonance gas with 76 mesonic and baryonic 
states up to masses of 2 GeV. The critical temperature, 
T c =170(180)MeV at SPS(RHIC) is roughly consistent 
with thermal-model fits to observed particle ratios (Hoj 
and predictions of lattice QCD [5l|. A freeze-out tem- 
perature of Tf ~ 120 MeV terminates the evolution and 
results in a total fireball lifetime of Tf o =10-12fm/c for 
central A-A collisions. The resulting temperature evolu- 
tion as a function of time r is displayed in Fig. [S] for SPS 
and RHIC. Note that there is little difference between 
mid (\y\ < 0.35) and forward rapidity (\y\ £ [1.2, 2.2]) for 
Au-Au collisions at RHIC due to the slowly varyin g ra - 
pidity density of charged particles over this y range [5 21 ] , 
cf. Ref. [53| for more details. 



9 



0.4 

0.35 

> 0.3 
O 

7J"0.25 
I 

b 0.2 

I 

H 0.15h 



0.05, 



i 


— RHIC-200 forward y: s/p B =155 : 

— RHIC-200 mid-y: s/p B =250 

.- SPS-17.3: s/p B =26 i 


: \ 


Central A-A : 


: > i 




= Oh 1 1 

rg ! MiX ! 









10 



Time (fm/c) 



FIG. 8: (Color online) Time profiles of temperature for central 
collisions of heavy nuclei (participant number A pa rt=380) at 
RHIC (\/s=200 AGeV; solid line: mid rapidity; dashed line: 
forward rapidity) and SPS (-^=17.3 AGeV; dot-dashed line). 



C. Rate Equation in Hot Medium 

We now proceed to the thermal rate equation to cal- 
culate the time dependence of the charmonium number, 
throughout the third stage of the heavy-ion reaction, 



dAT* 
~dT 



= -r$ iss (T) [N* - N e J(T)} 



(23) 



The rate equation is solved separately for 1 $>=J/ip, \c and 
if>', with initial conditions given by Eq. (fT!?|) . The loss 
term, — r^ ISS A^,, accounts for the dissociation of primor- 
dially produced charmonia and the gain term, T^N^ 1 , 
for the regeneration of charmonia via coalescence of c 
and c quarks. Both processes are governed by the 
temperature-dependent inelastic reaction rate, r^, lss (T), 
which is taken from Sec. Ill Bl (note that we do not em- 
ploy the p=0 value of the rate but rather a 3-momentum 
averaged value which is slightly larger; its precise value is 
obtained by matching the final yield of the loss term to 
the exact result obtained from solving the momentum- 
dependent Boltzmann equation (|31j) below). The equi- 
librium limit of each charmonium state, N^(T), is intri- 
cately related to their equilibrium properties, which will 
be discussed further below. To make the decomposition 
of the J/ij) number at any time r, 



(24) 



into (suppressed) primordial charmonia, N^ , and re- 
generated ones, A^ cs (r), more explicit, we exploit the 
linearity of the rate equation (f23|) . We define A^ rim (r) 
as the solution of the homogeneous rate equation, 



dA^" m 
dr 



-r$ iss Nf m 



(25) 



with the same initial condition as for the full rate equa- 
tion, N^ lm (0) = Ny(Q). The regeneration component, 
Ny S , then follows as the difference between the solution 
of the full and the homogeneous rate equation, which can 
be expressed as 



dA^ cs 
dr 



= -r 



diss 



(% eB - N C J) 



(26) 



with vanishing initial condition, N^ s (t < t* ) = 0. The 
onset time of regeneration processes, r*, is defined by 
T(t *) = T$ iss for each state 

Let us now return to the equilibrium limit of the char- 
monium abundances, N^(T), which we evaluate within 
the statistical model. Since the thermal production and 
annihilation rates of cc are believed to be small at SPS 
and RHIC energies, cc pairs are assumed to be exclu- 
sively produced in primordial N-N collisions and con- 
served thereafter. The open and hidden charm states are 
then populated in relative chemical equilibrium accord- 
ing to the canonical charm-conservation equation, 



1 N h(N op ) 

2 ° P Io(N op ) 



+ N hid 



(27) 



with N c5 : total number of charm-quark pairs from ini- 
tial production, N op = 7 c VpB n op : total number of all 
open-charm states with pertinent equilibrium density 
"op, Ahid = Tc^FB^hid- total number of all charmo- 
nium states with pertinent equilibrium density n^a, and 
7 C : charm-quark fugacity accounting for the deviation 
of chemical equilibrium with the heat bath (7 C =1 in 
full equilibrium). The ratio of modified Bessel func- 
tions, h (N op ) /J (N op ), on the right-hand-side of Eq. (127) 
is the characteristic canonical suppression factor which 
accounts for the exact conservation of net-charm num- 
ber, N c — N s , in the limit of small N c5 (canonical 
limit) (Mil: for N op < 1, one has Ii(N O p)/I (N O p) -> 
^N op , which acts as an additional (small) probability to 
enforce a vanishing net charm content in the system (i.e., 
both c and c have to be present simultaneously) . 

The open-charm number, A^p, is evaluated as fol- 
lows. For the QGP phase in the weak-binding scenario 
only charm quarks are counted as open-charm states. 
In the strong-binding scenario, the T-matrix calcula- 
tions of Ref. [IH suggest that cq and cq (charm-light) 
bound states survive in QGP up to ~1.3T C ; therefore, 
we count both charm quarks and the lowest-lying S- 
wave Z?-mesons (D, D* , D s and D*), as open charm 
states for T < 1.3T C . The charm-quark masses in the 
QGP correspond to the temperature-dependent ones dis- 
played in Fig. [TJ while for the meson resonances above 
T c we estimate from Ref. [ll| mo = raw — 2.0 GeV and 
mo 3 — tud* — 2.1 GeV (hyperfine splitting has been ne- 
glected) . For the HG phase all charmed hadrons listed by 
the particle data group (5(| are counted as open-charm 
states, with their vacuum masses. The hidden charm 
number, A^m, is evaluated in line with the existing char- 
monium states and their masses at given temperature T, 
but its contribution to N c5 is numerically negligible. 
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K(t) = 1 - exp(-r/r c oq ) (29) 
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FIG. 9: (Color online) Temperature dependence of the in- 
medium J/ip equilibrium limit using the statistical model 
in the QGP within the strong-binding scenario (dot-dashed 
lines: with and without D-meson resonances below and above 
1.3 T c ~234 MeV, respectively; solid line: smooth interpola- 
tion of the previous two cases; see text for details), the weak- 
binding scenario (dashed line) and in the HG for temperatures 
below T c =180MeV. 



Knowing n op , rihid and Vfb at each temperature, one 
can solve Eq. (pT7f for the charm-quark fugacity, 7 C (T), 
and apply it to compute the statistical equilibrium limit 
of each charmonium state as 



N$ = 7c ^FB 



(28) 



in terms of its equilibrium density, . In Fig. |H] we 
collect the numerical results of the statistical equilib- 
rium limit for J / ip abundances (excluding feeddown) for 
central 200AGeV Au-Au collisions at RHIC. The dis- 
continuity at 1.3 T c for the strong-binding scenario (dot- 
dashed line) is due to the inclusion of the D resonances 
in the QGP medium. We smoothly interpolate around 
the melting temperature for the Z?-mesons with a hyper- 
bolic tangent function (solid line) to represent a more 
gradual (dis) appearance of the D resonances (we have 
checked that this procedure has negligible impact on the 
calculation of observables in Sec. IIV|) . 

To achieve a more realistic implementation of the sta- 
tistical equilibrium limit, we apply two corrections to 
7V|, tat to schematically implement off-equilibrium effects 
of charm quarks in momentum and coordinate space. 
The former is aimed at simulating incomplete thermaliza- 
tion of the charm-quark p t spectra throughout the course 
of the thermally evolving bulk medium. It is expected 
that the coalescence rate from non- or partially thermal- 
ized c- and c-quar k sp ectra is smaller than for fully ther- 
malized ones [571 158| , since the former are harder than 
the latter and thus provide less phase-space overlap for 
charmonium bound-state formation. We implement this 
correction by multiplying the charmonium abundances 
from the statistical model with a schematic relaxation 



where r° q is a parameter which qualitatively represents 
the thermal relaxation time of charm quarks (it is one 
of our 2 main adjustable parameters in our phenomeno- 
logical applications in Sec. IIV[) . A rough estimate of this 
time scale may be obtained from microscopic calculations 
of this quantity within the same T-matrix approach as 
used here for charmonia, where the thermal char m- q uark 
relaxation time turns out to be r^ q ~3-10fm/c [l5l 159$. 
Such values allows for a fair description of open heavy- 
flavor suppression and elliptic flow at RHIC [HiJ E(| ■ The 
second correction is applied in coordinate space, based 
on the realization that, after their pointlike production 
in hard N-N collisions, the c and c quarks only have a 
limited time to diffuse throughout the fireball volume. 
At RHIC and especially at SPS only few cc pairs are 
produced (e.g., dN cS /dy ~ 1.2 in semicentral (6=7fm) 
Au-Au collisions at RHIC), and the hadronization time 
is smaller than the fireball radius. Thus, c and c will 
not be able to explore the full fireball volume but rather 
be restricted to a "correlation volume", V COTT [2l|, [sTJ] 
(the analogous concept has been successfully applied to 
strangeness production in p-A and A- A collisions in the 
SPS energy regime 62]). We implement this correction 
by replacing the fireball volume Vfb in the argument of 
the Bessel functions in Eq. ([2T)) by the correlation volume 
V COII j2l|, [24|. The latter is identified with the volume 
spanned by a receding cc pair, 



47T 

Korr(r) = —{r + (v c )t) 3 



(30) 



where tq ~ 1.2 fm represents an initial radius character- 
izing the range of strong interactions, and (v c ) is an aver- 
age speed with which the produced c and c quark recede 
from the production point; we estimate it from the aver- 
age pt in .D-meson spectra in p-A collisions [2l|, [6^, [64| 
as (v c ) ~ 0.55(0.6)c at SPS (RHIC). The correlation vol- 
ume leads to a significant increase of 7 C (since Iq/Ii is 
reduced) and thus of the modified ^ "equilibrium limit" 
due to an effectively larger cc density. 



D. Transverse-Momentum Dependence 

The transverse momentum dependence of in-medium 
charmonium production in A- A collisions has been sug- 
gested as a tool to better discriminate primordial and re- 
generated production [3(| ■ To investigate this observable 
within our approach we adopt the procedure suggested 
in Ref. [23j where the p t spectra for the 2 components are 
evaluated at hadronization based on the decomposition 
given by Eqs. ([24]) - ([26]) . For the primordial component 
we employ a more differential version of the rate equa- 
tion (|25j) . i.e., a Boltzmann transport equation, to de- 
scribe the evolution of the charmonium phase-space dis- 
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tribution functions, f^,(x,p, r), in a thermalized medium, 

p^d^f m (x,p,r) = r| iss (f,p,r) ff m (x,p,r) 

(31) 

where p a = = (m% + P 2 ) 1 ' 2 . The initial condition, 
fy Ilm (x,p,T ) = fy( x,T ) hip, t ), is still given as dis- 
cussed in Sec. IIII Al As mentioned before, we neglect 
elastic charmonium rescattering, and consequently the 
modification of the charmonium p t -spectra in the ther- 
mal stage is mostly due to the momentum dependence 
of the dissociation rates, T^ ss (x,p, r), given in Eq. (fTO]) 
and plotted in the lower panel of Fig. [3J We also account 
for "leakage" effects, i.e., charmonia escaping the fire- 
ball volume are no longer subject to suppression, which 
leads to the opposite trend compared to the increased 
suppression with p induced by the momentum depen- 
dence of the dissociation rate (charmonia with low pt are 
more strongly suppressed since they stay longer within 
the fireball). At the moment of freeze-out, T{ , we obtain 
the final p t spectra of the primordial component by inte- 
grating the spatial part of the solution of the Boltzmann 
equation, 



dm 



Ptdpt 



= 2n / d xfy(x,p t ,T fo ) . 



(32) 



A microscopic evaluation of the momentum depen- 
dence of the regeneration component using the quasifree 
mechanism requires to compute a 3-to-2 process, i + c + 
c — > i + "J. In the framework of the Boltzmann equa- 
tion, this is a rather involved calculations which will be 
reported in an upcoming publication [65]. For the time 
being, we approximate the pt distributions of regenerated 
charmonia by local thermal distributions boosted by the 
transverse flow of the medium, amounting to a standard 
blastwave description (66j . 



Ptdpt 



oc m t 



rdrKi 



mt cosh y t 
T 



la 



p t sinh y t 
T 



(33) 

(mt — \frn%, + p 2 )■ The medium is characterized by 
the transverse-flow rapidity yt = tanh _1 u t (r) using a 



linear flow profile v t (r) 



with a surface velocity 



v s = diT m i X and transverse fireball radius R — i?(r m i x ) 
as given by the fireball expansion formula, Eq. (|2"2"j) . at 
the end of the mixed phase, r m j x . We evaluate the blast- 
wave expression at the hadronization transition (T c ) and 
neglect rescattering of VP's in the hadronic phase. Two 
additional effects are neglected in this treatment which 
to a certain extent tend to compensate each other: on 
the one hand, due to incomplete charm-quark thermal- 
ization, one expects the regenerated charmonium spectra 
to be harder than in equilibrium, but, on the other hand, 
a good part of the regeneration occurs before the mixed 
phase [H,[63] (see also Fig. ([14)) below) so that the evalu- 
ation of the blastwave expression at the end of the mixed 
phase presumably overestimates the blue shift due to the 
flow field. An explicit evaluation of the gain term with 



realistic (time-dependent) charm-quark spectra within a 
Boltzmann transport equation [65] will be able to lift 
these approximations. 

The total charmonium p t spectra are the sum of the 
suppressed primordial and regenerated parts, Eqs. (I32[) 
and (|33p . with an absolute normalization following from 
the decomposition of the pt-integrated rate equation, 
Eq. (|2"4"1) . at thermal freezeout. The average p t squared 
is readily computed as 



(Pt ) prim,rcs — 



I d 2 p t P 2 



Ptdpt 



Jd 2 Pt 



(34) 



pt dpt 



which we compare to experimental data as part of the 
following section. 



IV. COMPARISON TO SPS AND RHIC DATA 

In this section we present and discuss the numeri- 
cal applications of the above framework to J/tp data in 
URHICs at SPS and RHIC. For each observable, we con- 
front the results of the strong- and weak-binding scenario 
in an attempt to discriminate qualitative features. The 
feeddown to J/ip from \c and ip' states is taken into ac- 
count, assuming fractions of 32% and 8%, respectively, 
for primordial production in pp collisions. We have di- 
vided the discussion into the centrality dependence of in- 
clusive J/tp yields in Sec. IIV Al and their p t dependence 
in Sec. IWBl 



A. Centrality Dependence of Inclusive Yields 

The J I tp yield in A- A collisions is usually quantified in 
terms of the nuclear modification factor as a function of 
centrality, 



Raa(P) 



(35) 



where N co n(b) is the number of binary collisions of the 
incoming nucleons at impact parameter b. Before we 
turn to the results, we recall the two main parameters 
in our approach, which are the strong coupling constant, 
a s , the thermal charm-quark relaxation time, r® q The 
former controls the inelastic charmonium reaction rate 
and the latter the magnitude of the $ equilibrium lim- 
its. We adjust them to approximately reproduce the in- 
clusive J/tp yield for central A- A collisions at SPS and 
RHIC, within reasonable bounds. For a s we find that a 
common value of 0.32, which is at the upper end of the 
value in the Coulomb term in the QQ free energy, can be 
used, in combination with T^r q =3.8fm/c for the strong- 
binding scenario and T° q =1.6fm/c for the weak-binding 
scenario. For simplicity, we refrain from introducing an 
additional temperature dependence into these parame- 
ters. The composition of the total yield, its centrality 
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FIG. 10: (Color online) Results of the thermal rate-equation 
approach for J/ip production (normalized to Drell-Yan pairs) 
versus centrality at SPS, compared to NA50 data [6a . [69l |. 
Solid lines: total J/ip yield; dashed lines: suppressed primor- 
dial production; dot-dashed lines: regeneration component; 
dotted lines: primordial production with CNM effects only. 
Upper panel: strong-binding scenario; lower panel: weak- 
binding scenario. 



FIG. 11: (Color online) Results of the thermal rate-equation 
approach for the nuclear modification factor versus centrality 
at mid rapidity at RHIC, compared to PHENIX data [70l ]. 
Solid lines: total J/ip yield; dashed lines: suppressed primor- 
dial production; dot-dashed lines: regeneration component; 
dotted lines: primordial production with CNM effects only. 
Upper panel: strong-binding scenario; lower panel: weak- 
binding scenario. 



dependence and pt spectra can then be considered as a 
prediction within each of the 2 scenarios. 

We begin with J/ip production in ^=17.3 AGeV Pb- 
Pb collisions at SPS, for which we compare our results in 
the strong- and weak-binding scenario with NA50 data in 
Fig. [TO] For these data, the denominator in Eq. (f35]) is 
replaced by the number of Drell-Yan dileptons at high 
mass, while the numerator includes the branching ra- 
tio into dimuons. The pertinent proportionality factor, 
equivalent to the pp limit of this ratio (47. Oil. 4 [71j), and 
the CNM-induced suppression (dotted line in Fig. [T0|) are 
inferred from latest NA60 p- A measurements 43] , which 
we reproduce using the Glauber model formula, Eq. (fT5|) , 
with <7 a bs=7.3mb. The suppression of the primordial 
component (dashed line) relative to the nuclear absorp- 
tion (dotted line) represents the "anomalous" suppres- 
sion by the hot medium, which increases with centrality 



due to higher initial temperatures and longer fireball life- 
times. The regeneration component increases with cen- 
trality as well, mostly due to the increase of the IZ-i actor 
and the larger lifetime which facilitates the approach to 
the equilibrium limit according to Eq. (|26p. According 
to detailed balance between dissociation and regenera- 
tion, an increase in the former also implies an increase 
in the latter. The sum (solid line) of primordial and re- 
generation contributions describes the centrality depen- 
dence of the inclusive J/ip yield at SPS reasonably well 
in both scenarios. In the strong-binding scenario the pri- 
mordial component is dominant and the majority of the 
anomalous suppression originates from the dissociation 
of Xc and ip', since at the temperatures realized at SPS 
(T ~ 200 MeV) the quasifree dissociation rates for \c 
and ip' are much larger than those for J/ip, recall Fig. [3] 
In the weak-binding scenario, however, the regeneration 
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FIG. 12: (Color online) Results of the thermal rate-equation 
approach for the nuclear modification factor vs. centrality at 
forward rapidity compared to PHENIX data [701 ] . Solid lines: 
total J/ip yield; dashed lines: suppressed primordial produc- 
tion; dot-dashed lines: regeneration component; dotted lines: 
primordial production with CNM effects only. Upper panel: 
strong-binding scenario; lower panel: weak-binding scenario. 



yield becomes comparable to the primordial one for semi- 
/ central collisions due to larger dissociation rates and the 
smaller charm-quark equilibration time scale. 

Next we examine the centrality dependence of J/ip pro- 
duction in 200AGeV Au-Au collisions at RHIC, first fo- 
cusing on mid rapidity (\y\ < 0.35), as shown in Fig. 1111 
The suppression due to CNM effects (dotted line in 
Fig. [TTj) is inferred from latest PHENIX d-Au measure- 
ments, which we reproduce using the Glauber model for- 
mula Eq. JISJ with cr abs =3.5mb. For N part ~ - 100, 
the composition of primordial and regeneration contribu- 
tions is quite comparable to the SPS for Ar part ~ — 400 
within both scenarios. Beyond N paTt ~ 100, suppression 
and regeneration continue to increase, leveling off at an 
approximately 50-50% (20-80%) partition for primordial 
and regeneration in the strong-binding (weak-binding) 
scenario in central collisions. 

Let us now turn to J / ip production at forward rapidity 



FIG. 13: (Color online) Results of the thermal rate-equation 
approach for the ratio of Raa for J/ip at forward and mid 
rapidity versus N par t in strong (solid line) and weak (dashed 
line) binding scenarios compared to PHENIX data [7(j|. In 
the upper two curves, CNM effects have been divided out in 
both numerator and denominator of the ratio. 



(\y\ € [1.2, 2.2]) at RHIC, shown in Fig.[TJ Again, both 
strong- and weak-binding scenarios reproduce the exper- 
imental data fairly well, with similar relative partitions 
for primordial and regeneration contributions as at mid 
rapidity. However, one of the "puzzles" about J/ip pro- 
duction at RHIC is the fact that the total J/ip yield is 
more strongly suppressed at forward rapidity than at mid 
rapidity. In our approach this follows from the stronger 
shadowing at forward rapidity leading to less primordial 
production for both J/ip and cc pairs. The former (latter) 
leads to a reduction of the direct (regeneration) compo- 
nent. Since the thermodynamic properties of the fire- 
ball are quite similar at mid and forward rapidity (recall 
Fig. , charmonium suppression and regeneration in the 
hot medium are very similar between the mid and for- 
ward rapidity as discussed in Ref. [Hj]. To quantify the 
difference at forward rapidity and mid rapidity we display 
the ratio between the corresponding Raa's in Fig. 1131 for 
both scenarios, which clearly illustrates the importance 
of CNM effects to properly reproduce the data. 

In Sec. Ill Al we have argued that the strong- and weak- 
binding scenarios discussed here may be considered as 
limiting cases for J/ip binding in the QGP, as bracketed 
by the identification of the heavy-quark internal and free 
energies with a Q-Q potential. From the results above 
we believe that these scenarios also provide a reasonably 
model- independent bracket on the role of suppression and 
regeneration effects, in the following sense: At SPS, the 
strong-binding scenario defines a "minimal" amount of 
dissociation required to provide the anomalous suppres- 
sion beyond CNM effects (a small regeneration compo- 
nent is inevitable due to detailed balance). The applica- 
tion to RHIC energy then implies an approximately equal 
partition of primordial and regenerated charmonia in cen- 
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mixed phase, since the HG reaction rates are small. In 
the weak-binding scenario the time-dependent J/ip yield 
exhibits a "dip" structure around r ~ 1.5fm/c because 
the large dissociation rates suppress primordial J/ ip very 
rapidly and regeneration only starts after the medium 
temperature falls below the J/ip dissociation tempera- 
ture ft — 1.25 T c ). This scenario thus is closest in 
spirit to the statistical hadronization model [l9( where 
all initial charmonia are suppressed (or never form to be- 
gin with, except for corona effects) and are then produced 
at the hadronization transition. 
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FIG. 14: (Color online) J/ip abundance as a function of time 
in central (Af part =380) Au-Au collisions at RHIC. Solid lines: 
total; dashed line: primordial component; dot-dashed line: 
regeneration component; double dot-dashed line: equilibrium 
limit of J/ip, JVju. Upper panel: strong-binding scenario; 
lower panel: weak-binding scenario. 



tral Au-Au, not unlike Ref. [22J where the vacuum char- 
monium binding energies ( "strong binding" ) have been 
used in the QGP (together with the gluo-dissociation) . 
On the other hand, in the weak-binding scenario, a large 
part of the J/ip yield in central A- A is due to regen- 
eration even at SPS, limited by the constraint that for 
sufficiently peripheral collisions (and at sufficiently large 
pt) a transition to primordial production compatible with 
p-A data should be restored. Clearly, for central A- A at 
RHIC (and certainly at LHC) the final yield is then dom- 
inated by regeneration. Since both scenarios describe the 
inclusive yields reasonably well, it is mandatory to inves- 
tigate more differential observables to find discriminating 
evidence which will be pursued in the following section. 

It is instructive to examine the time evolution of J/tp 
production in the two scenarios, displayed in Fig. [14] 
for central collisions at mid rapidity at RHIC (excluding 
feeddown from \c and ip 1 ). In both scenarios most of the 
dissociation and regeneration occur in the QGP and the 



B. Transverse-Momentum Spectra 

The results of the previous section suggest that, within 
the current theoretical (e.g., charm-quark relaxation time 
T° q ) and experimental uncertainties both of the "limit- 
ing" scenarios can reproduce the centrality dependence 

of the inclusive R J J^ '(N part ) reasonably well at both SPS 
and RHIC energies. However, the composition between 
suppression and regeneration yields is rather different 
which ought to provide a key to distinguish the two sce- 
narios. The obvious "lever arm" are charmoniump t spec- 
tra [301 ] . One expects that the primordial component is 
characterized by harder p t spectra (following a power law 
at high p t ) while the regeneration component produces 
softer p t spectra characterized by phase-space overlap of 
(partially) thermalized charm-quark spectra. However, 
in practice, the transition from the "soft" recombination 
regime to the "hard" primordial regime is quite uncer- 
tain; e.g., collective flow and incomplete thermalization 
of c-quarks can lead to a significant hardening of the re- 
generated J/ip spectra, while a dissociation rate which 
increases with 3-momentum [231 ] can induce a softening 
of the spectra of the surviving primordial charmonia. 

For a concise discussion of the pt dependence of J/tp 
as a function of centrality at SPS and RHIC we here fo- 
cus on the average p\ , as compiled in Fig. fTS] At the 
SPS (left panels), the centrality dependence of {p\) is 
largely dictated by the the Cronin effect in the primor- 
dial component, especially in the strong-binding scenario 
where this contribution dominates the yield at all cen- 
tralities. The momentum dependence of the dissociation 
rate induces a slight suppression of (pf) at large central- 
ity compared to the case where only CNM effects are in- 
cluded (dashed vs. dotted line) (2^|. In the weak-binding 
scenario, larger contributions from regeneration induce a 
slight "dip" structure at intermediate centralities due to 
a rather small collective flow at the end of the (relatively 
short) mixed phase in these collisions. 

At RHIC energy (middle and right panels of Fig. fTS]) , 
the individual primordial and regeneration components 
show qualitatively similar behavior for (pf){N pav t) as at 
SPS, i.e., an increase due to Cronin effect and collective 
flow, respectively. At mid rapidity, the general trend 
is that with increasing centrality the growing regenera- 
tion contribution pulls down the average (pt), in qual- 
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FIG. 15: (Color online) Results of the thermal rate-equation approach for (pi ) vs. centrality at SPS (left panels, compared 
to NA50 data [38l. l72l]) and RHIC at mid and forward rapidity (middle and right panels, respectively, compared to PHENIX 
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primordial component (dashed line), the regeneration component (dash-dotted line) and primordial production with CNM 
effects only (dotted lines). The upper (lower) panels correspond to the strong- binding (weak- binding) scenario. 



itative agreement with the data. The curvature of the 
(pi) (iV part ) dependence, which appears to be negative in 
the data, is not well reproduced, neither by the strong- 
nor by the weak-binding scenario, even though the de- 
viations are smaller in the former. A more microscopic 
calculation of the gain term, together with more accurate 
estimates of the Cronin effect, are warranted to enable 
more definite conclusions. For both rapidity regions, the 
(pi) of the suppressed primordial component is slightly 
larger in the weak- than in the strong-binding scenario. 
This is caused by the stronger 3-momentum increase of 
the dissociation rate in the strong-binding scenario, recall 
lower panel of Fig. [3] 

The overall comparison to SPS and RHIC data for the 
pt dependence of J/ip's seems to indicate a slight prefer- 
ence for the strong-binding scenario. This is mostly de- 
rived from the observation that for peripheral collisions 
the experimentally observed (pi) essentially follows the 
extrapolation of the Cronin effect, suggesting J/ip pro- 
duction of predominantly primordial origin (the collective 
flow imparted on the regeneration component appears to 
be too small at these centralities). 

The average pi mostly characterizes the momentum 
dependence of charmonium production at low and mod- 
erate pt where most of the yield is concentrated. Recent 
RHIC data [z3>[2H have triggered considerable interest in 
J / ip production at high p t ~ 5—10 GeV which is expected 
to provide complementary information. It was found 

that the suppression in R'^(p t > 5 GeV) in Cu-Cu col- 
lisions is reduced compared to the \ow-p t region, with 



.R^A-values of —0.7-1 or even larger (74|. This is quite 
surprising in light of the light-hadron spectra measured 
thus far at RHIC which all exhibit stronger suppression of 
Raa — 0.25 for p t > 6 GeV (even electron spectra from 
open heavy flavor, i.e., charm and bottom decays). It 
also appears to be at variance with the thermal J/ip dis- 
sociation rates which, if anything, increase with momen- 
tum (recall Fig. [3]) and thus imply a stronger suppression 
at higher p t . Furthermore, the leakage effect referred to 
above is not strong enough to produce the experimentally 
observed increase in 

r aa(pt > 5 GeV) However, 
as has been pointed out more than 20 years ago [75ll76j . 
the finite formation time of the charmonium states, in 
connection with Lorentz time dilation at high momen- 
tum, can lead to a reduction in the effective absorption 
cross section (or thermal rate) as long as the bound-state 
wave function has not developed (the pre-resonance state 
is more compact than the fully formed J/ip). Further- 
more, for high-pi J/ip spectra significant feeddown con- 
tributions from the decay of B-mesons [73] is expected 
whose total yield is not suppressed. In Ref. [78| we have 
schematically implemented these effects into our rate- 
equation approach (as in the present paper) and found 

that i?^ -values of up to -0.85 at p t ~ 10 GeV [z| 
can be recovered in Cu-Cu. Here, we extend these cal- 
culations to provide predictions for central Au-Au colli- 
sion and specifically address the question whether high-p t 
J/ip production can help to disentangle the strong- and 
weak-binding scenarios. Here we do not readjust any pa- 
rameters relative to the above calculation which implies 
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FIG. 16: Nuclear modification factor for J/ip pt spectra 
in central 200 AGeV Au-Au collisions including formation- 
time effects (fte) and B-meson feeddown (Bfd) contributions. 
PHENIX data [7(j] are compared to our rate-equation calcu- 
lations in the strong- and weak-binding scenarios (upper and 
lower panel, respectively). 



an increase of 15-20% in the total yield. Our results for 
R J j[^{pt) are displayed in Fig.flrJlup to p t = 10 GeV. We 
find that the suppression is reduced to about 0.5 at the 
highest pt, compared to about 0.4 at low p t . This is sim- 
ilar to the moderate enhancement found in the Cu-Cu 
case (where it went from from ^0.65 to 0.85). More sur- 
prisingly, the high-p t suppression is very similar in both 
strong- and weak-binding scenarios, despite the fact that 
the high-pt yield is exclusively due to the primordial com- 
ponent whose strength is very different in the 2 scenarios 
at low pt- The reason is the 3- momentum dependence 
of the dissociation rates, which become quite similar in 
the 2 scenarios at large 3-momentum: at p ~ 10 GeV, 
the difference in the energy-threshold due to binding en- 
ergies of several 100 Me V becomes less relevant so that 
a collision with almost any thermal parton is energetic 
enough for dissociating the bound state. 

Let us finally comment on the elliptic flow, V2(j>t), of 



the J/ip, which was hoped to be another good discrim- 
inator of primordial and regenerated production. For 
the former, a nonzero i>2 is basically due to the path- 
length difference when traversing the azimuthally asym- 
metric fireball, typically not exceeding 2-3%. For the 
latter, much larger values can be obtained if the coalesc- 
ing charm quarks are close to thermalized. However, as 
pointed out in Ref . [78] , the well-known mass effect sup- 
presses the V2(pt) for heavy particles at pt < m; it is 
precisely in this momentum regime where the regenera- 
tion component is prominent. Thus, we predict that in 
both strong- and weak-binding scenarios the total J/tp 
V2(j>t) does not exceed ~3% at anyp t . The only alterna- 
tive option we can envision are strong elastic interactions 
of the J/ip which is only conceivable in the strong bind- 
ing scenario to avoid break-up in scattering off thermal 
partons [64j |. 



V. CONCLUSION 

The main goal of this work was to construct charmo- 
nium spectral functions which are constrained by thermal 
lattice-QCD computations, and to apply them to experi- 
mental data in heavy-ion collisions. Employing a thermal 
rate equation we have implemented equilibrium proper- 
ties of charmonia (binding energy, constituent charm- 
quark mass and dissociation rate) as extracted from 
a thermodynamic T-matrix calculation. The resulting 
spectral functions for two "limiting" scenarios, with small 
and large J/ip dissociation temperatures, have been ver- 
ified to produce a weak temperature dependence in per- 
tinent euclidean correlators, roughly compatible with 
1QCD. We have argued that these two scenarios may 
serve as generic representatives for charmonium kinet- 
ics in heavy-ion collisions, bracketing strong and weak 
in-medium binding. Therefore, we believe that the qual- 
itative conclusions drawn from these two scenarios should 
be rather model-independent. Within current theoretical 
uncertainties (especially for the degree of charm-quark 
relaxation and its impact on the J / ip regeneration yield) 
and essentially two fit parameters (r c and a s controlling 
the regeneration yield and suppression strength) both 
scenarios can reproduce SPS and RHIC data for the cen- 
trality dependence of inclusive J/ip production reason- 
ably well. However, the partition of primordial and re- 
generated yields is quite different in the two scenarios: 
the former dominates for strong binding (down to 50% 
in central Au-Au at RHIC), while for weak binding re- 
generation largely prevails at RHIC energies (except for 
peripheral collisions). We have investigated the p t depen- 
dence of the J/ip yield in both scenarios and found that 
differences in the average p\ reach up to 20% in semicen- 
tral Au-Au collisions. The strong-binding scenario seems 
slightly favored in this observable, but theoretical uncer- 
tainties (e.g., the blast- wave treatment of the regenera- 
tion component) prevent us from definite conclusions at 
this point. We also have to await more accurate A- A 
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and p(d)-A data, where the latter determine the input 
in terms of shadowing, nuclear absorption and Cronin 
effects. 

Further developments of the theoretical approach are 
in order, and we plan to improve our calculations in 
several respects. First, as mentioned before, work is in 
progress employing a Boltzmann transport equation (|31[) 
extended by a microscopic treatment of the gain term 
to treat charmonia regeneration and dissociation on the 
same footing. This will enable an explicit account of 
(time-dependent) charm-quark phase-space distribution 
in c-c recombination reactions, as following, e.g., from 
realistic Langevin simulations [59[ with constraints from 
the T-matrix formalism and from open-charm observ- 
ables. Second, hydrodynamic simulation of the medium 
evolution could be employed for a more detailed and re- 
alistic description of the temperature and the flow field 
of the underlying medium, especially in coordinate space. 
Third, a microscopic model for primordial cc and char- 
monium production is warranted to better disentangle 
nuclear shadowing and absorption in the pre-equilibrium 
stage, including formation-time effects. This would im- 
prove the initial conditions to the kinetic approach in 



the hot medium. These developments will eventually 
produce a comprehensive approach which can serve as a 
quantitative bridge between theoretical studies of char- 
monia in the QGP (and HG) and their phenomenology in 
heavy-ion collisions. As our investigations in the present 
paper have shown, quantitative studies at the 10-20% 
level are needed to deduce basic properties of the strong 
force (such as color screening of Coulomb and confining 
interactions) from heavy-ion collisions. These insights 
are also pivotal to improve our knowledge of the phase 
structure of hot and dense QCD matter. 
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